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Abstract 

We present a formulation of the shallow water equations that em- 
phasizes the conservation of potential vorticity. A locally conservative 
semi-Lagrangian time-stepping scheme is developed, which leads to a 
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system of three coupled PDE’s to be solved at each time level. We 
describe a smoothing analysis of these equations, on which an effective 
multigrid solver is constructed. Some results from applying this solver 
to the static version of these equations are presented. 

1 Formulation of the Shallow Water Equa- 
tions 


The shallow water equations provide a two-dimensional prototype of the 
equations needed for three-dimensional simulations of atmospheric motions 
[l] [2]. They are useful for testing the viability of new numerical schemes 
for atmospheric simulation because they share many of the properties with, 
but lack the full complexity of, a full three-dimensional system. The shallow 
water equations can be written as 


du 

dt 


-4>X + fv, 


dv 

dt 


= -<t>y ~ /«, 


dt 


-4>D, 


( 1 ) 

( 2 ) 

( 3 ) 


where u and v are the velocity components of the wind, D = u x + v y is 
the divergence of the velocity, / is the Coriolis parameter, and <f> is the 
geopotential height, assumed to be a positive function. The derivatives are 
material derivatives, that is, 


d d d d 
Jt ~ u lh + v fy + di' 


( 4 ) 


A considerable amount of effort has gone into designing numerical methods 
that will solve these equations (see for example the references cited in [1]). 
The purpose of this paper is to study a multigrid scheme applied to a form 
of these equations that is of special physical interest. 

There are many possible formulations of the shallow water equations. 
We will derive a different formulation from the one above that has certain 
physical and numerical advantages. To this end, we define vorticity by 


C = V X - Uy. 


( 5 ) 
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Then, subtracting the ^-derivative of Eq. 2 from the ar-derivative of Eq. 1 
gives 

J t « + /) = -« + f)D. ( 6 ) 

Solving for D in Eq. 3 and substituting it into Eq. 6 yields 

Eq. 7 is important in practice because it clearly asserts that the physical 
quantity (( + /)/<?!>, called potential vorticity, is conserved in time along any 
Lagrangian trajectory. 

Now adding the x-derivative of Eq. 1 and the y-derivative of Eq. 2 gives 

i£ = -VV-V.(/kxV)-W, (8) 

at 

where k = (0,0, 1), V = (u,u,0), and N = (u x ) 2 + (uy) 2 + 2 v x u y . It is not 
hard to see that Eqs. 3, 7, and 8 are equivalent to the original formulation 
of the shallow water equations (Eqs. 1-3), but they are not yet in the form 
we wish to consider. 

From the point of view of a multigrid solver, we will see that it is conve- 
nient to rewrite these equations in terms of the geopotential, <j), the stream 
function, t/’, and the velocity potential, x • The latter two variables satisfy 

V = k x + Vy, (9) 


C = vv, (10) 

D = V 2 X . (11) 


Using these variables, we arrive at the form of the shallow water equations 
used in this paper: 


£[Y!i±Zj = o. 

dt 1 <f, 1 


( 12 ) 


dV^x 

dt 


—V 2 <f) — V.(/k x Vx — /V^) — N, 


(13) 


d^_ 

dt 


-<^V 2 x, 


(14) 
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where N — (0 xy Xxx) H“ ( 0xy H~ Xyy) 2(0 yy Xxy){^xx ~f" Xxy)* 

These equations have several attractive properties. As already noted, they 
emphasize the conservation of potential vorticity along Lagrangian trajecto- 
ries. Furthermore, we shall show in Section 2 that, when a semi- Lagrangian 
approach is taken for the time derivatives, all of the variables appear in po- 
tential form in the resulting equations. This means that a simple vertex 
centered grid is sufficient to discretize the problem spatially; a staggered grid 
is not needed. This fact should be particularly useful when the problem is 
posed on a spherical domain. Finally, we shall see in Section 3 that these 
equations are well suited for multigrid solution. 

An ideal domain for simulating atmospheric motions is a sphere. However, 
a spherical coordinate system introduces many difficulties that may confuse 
the task of developing an efficient solver for the equations at hand. Thus, as 
a first step in determining the feasibility of applying multigrid methods to 
our formulation of the shallow water equations, we have chosen to solve the 
system on a cylindrical domain. Specifically, we consider a domain that is 
periodic in the x direction with length d and includes y in the range [0, L]. 
We set x = 0 = 0 and <f> = <f> 0 at the y boundary, where <f > o is a given 
constant. We assume that the Coriolis parameter may be written as 

f = fo + Py, (15) 

with f 0 and (3 constants. This model allows us to determine the effective- 
ness of multigrid methods for these equations without the complications of 
constructing a full three-dimensional global atmospheric model. 

2 A Semi-Lagrangian Time Stepping Scheme 

Eqs. 12-14 are written in a Lagrangian reference frame in which the evo- 
lution of the fluid is observed along the paths of imaginary fluid particles. 
There are some obvious disadvantages of evolving a set of particles along 
Lagrangian trajectories numerically. In particular, a grid that is initially 
uniform will in general become very irregular, often leading to a degradation 
of global accuracy. As a compromise, semi-Lagrangian methods have been 
developed to produce numerical methods that preserve the advantages of 
regular grids while simultaneously taking advantage of the Lagrangian form 
of the equations. There is an extensive body of literature describing these 
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Figure 1: A schematic diagram showing the main quantities used in the 
calculation of the departure points for the semi-Lagrangian time-stepping 
scheme. The exact trajectory is represented by a solid line and the approxi- 
mate trajectory with a dashed line. 

methods. In particular, [1] provides an excellent review of the application 
of semi-Lagrangian methods to meterological problems. This reference de- 
scribes in detail a semi-Lagrangian scheme for the integration of Eqs. 1-3. 
The scheme we describe in this section is an adaptation of this scheme to 
our reformulation of the shallow water equations, and the reader is urged to 
consult [1] for more detail. 

The fundamental idea of a semi-Lagrangian scheme is to impose a regular 
grid at the new time level, and to backtrack the fluid trajectories to the 
previous time level. At the old time level, the quantities that are needed 
are evaluated by interpolation from their known values on a regular grid. In 
general, as is the case in our problem, the velocity field at the new time step 
is unknown, so the critical problem in this idea is the computation of the 
trajectory departure points. 

A schematic representation of the quantities involved in computing the 



departure points is shown in Fig. 1. The displacement between a grid point 
on the new time level, x m (t), and the departure point of the trajectory leading 
to this point on the previous time level, x^,(f — Af), is denoted by a m . If the 
velocity field is considered to be constant from t — At to t , then a m satisfies 
the equation 

3 A/ 

a m = A tV(x m - - — ). (16) 

The velocity at time t — Ai/2 may be defined by extrapolation from the two 
previous time levels by 

V(x, < - y) = ^V(x. ( - At) - iv(x, t - 2A() + C>(Ai 2 ). (17) 

Eqs. 16 and 17 give an implicit equation for a m in terms of the known velocity 
field at two previous time levels, and we may consider an iterative method 
for determining the correct a m . Assuming that a suitable approximation is 
made, then x m — a m /2 would not generally lie on a grid point, so the velocities 
at this point must be obtained by interplation. It has been shown [4] [5] [6] 
that for problems of this type it is sufficient to use linear interpolation to 
define the quantities in Eq. 17. It is also known [7] that succesive iteration 
for the solution of Eq. 16 converges provided 


At < 


1 

max(|ti„i, |u,|, |»,|, |d,[ ' 


(18) 


Once the a m are known, the departure point values of the variables in 
our equations are defined as illustrated by 


- At) = <f>(x m - a m ,f - At). (19) 

Again, these values must be interpolated from known values at the grid 
points. It has been found [4] [5] [6] that it is advantageous to do this using 
cubic interpolation. A material time derivative may then be discretized by 

Tt = h [m - f(t ~ Ai)l> (20) 

and nonderivative quantities can be represented by the simple average 

* = \w) + nt - At)]- (2i) 
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Using this discretization, our formulation of the shallow water equations 
may be manipulated to show that the equations that determine the solution 
of the system at a new time level are 

Vyif-//=0, (22) 


W + r[W + + Pxt - Wt - /VV + ] = h, 
^ + [1 + rV 2 x + ] = /3, 


(23) 

(24) 


where 


, v V + /* 

* l 4>* 

/ 2 = V.[V*-r(/kxV' + Vf)], 

h = ^(i - rVV), 


(25) 

(26) 
(27) 


and r = Af/2. The starred quantities are evaluated at the trajectory de- 
parture points at the previous time level, and the superscript -f refers to 
quantities defined on a regular spatial grid at the new time level. We refer 
to Eqs. 22-24 as the static equations. The superscript + will be omitted in 
what follows. 

The numerical algorithm needed to integrate our form of the shallow 
water equations splits naturally into two pieces. The first task is to compute 
the departure point quantities needed to define /i, fo, and fo. This is done 
in the manner outlined above, using information from two previous time 
levels. The velocity field at any time level may be obtained from x and 0 
using u = -i/f y + Xx and v = rf> x + Xy 0 nce the departure point quantities 
are known, the second task is to solve the static equations. As we shall 
demonstrate below, it is possible to construct an efficient multigrid solver for 
these equations. Note that nowhere in this method is it necessary to solve 
Eqs. 10 and 11 for x and xp in terms of u and v. 


3 Coupling Analysis of the Static Equations 

The coupling between the equations in any system of equations plays a piv- 
otal role in the behavior of the system. In particular, when discretized sys- 
tems of PDE’s are to be solved by multigrid, the coupling of the equations 



must determine the character of the relaxation schemes that are to be ap- 
plied. Fortunately, a straightforward method for analyzing the coupling of 
a system and its relation to constructing a multigrid solver is available [8]. 
In this section we apply this method to the static equations derived above. 
Throughout the section, we use the definitions and notation of [8], to which 
the reader is referred for an understanding of the technique we are about to 
use. 

The linearized static equations are given in brief as follows: 

( -/rW - T0d y rV‘ V 2 +°r fid, ) [ \ ) = ( /). (28) 

V 0 1 + I \ x } l h ) 

In constructing this system, we have associated variables with equations in 
the natural way; that is, i/>, <f>, and x are associated with Eqs. 22, 23, and 
24, respectively. 

The order array and weight array for this system are 

' 2 0 N ' 

Q= 2 2 2 , (29) 

TV 0 2 

and 

' N 2 N ' 

W = 0 JV 0 , (30) 

_iV 2 1V_ 

respectively. 

To account for finite mesh size effects, we need the scaled coefficient array 


'1 -h N ' 

C= -/ 1 r_1 (31) 

N 1 

L T<P J 

The computation of these arrays is straightforward. The method of [8] is 
almost automatic, and the arrays are included here explicitly only for com- 
pleteness. From these arrays, the coupling graph may be constructed, as 
shown in Fig. 2. 
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Figure 2: The coupling graph for the static equations. The finite mesh- 
size coupling coefficients are /i (1 -+ 2), / (2 -» 1), 1/ T ( )> an 

[1 + rV 2 x]/r<£ (3 - 2). 

We may conclude immediately from the coupling graph that Eq. 2 
weakly coupled to Eq. 23 when 

/,/ft 2 « i, ( 32) 

and that Eqs. 23 and 24 are weakly coupled when 

( 1 + tV 2 x)^ <<1 - 

This implies that it both of these conditions are satisfied, then each equation 
may be relaxed separately, as though the system were ful y ec ° u P e 

We now need to estimate the quantities in these coupling conditions using 
a physically realistic solution of a slightly different version of the shallow 
water equations. The equations we are dealing with assume that the surface 
of the fluid is free. To fix the surface profile of the fluid (the so-called ngi 
lid’ condition), we set d<f>/dt = 0 in Eq. 1. It can then be shown by direct 
substitution that the following is an exact form for the resulting Rossby- 
Haurwitz wave solution: 

(34) 


u 


= U — Al cos ly s\n k(x — ct), 



f 

fo = 1 X 10 -*s~ 1 ,jl = 1.57 x 10 -11 m -1 5 _1 

r 

5005 

d 

1 x 10 7 m 

h 

~ 10~ 8 

L 

5 x 10 6 m 

y*2 

~ 10- 7 

<f> 0 

1 x 10 4 m 

h 

~ 10 4 


Table 1. Some typical physical parameters for the shallow water equations. 

v = Ak sin ly cos k(x — ct), ( 35 ) 

$ ~ do — foUy — ~/ 3 Uy 2 + A sm k(x — ct)[f sin ly — (c — U)l cos ly] + 

1 


jA 2 [r cos 2 k(x — ct) + k 2 cos 2/y], 


(36) 


where A and U are constants, c — U — (/3 2 /(k 2 + Z 2 )) is the Rossby-Haurwitz 
phase speed, k = < l'Kvnj L for integer m, and l = mr /d for integer n , Waves of 
this type are the dominant feature of large scale weather motions. This solu- 
tion satisfies different boundary conditions from the problem we are treating, 
but it is nevertheless useful for estimating the size of the parameters in our 
system. It can be shown for this solution that 


V 2 x = 0 (37) 

and 

V 2 0 = -A(k 2 + l 2 ) sin ly sin k(x - ct). (38) 

Some typical numerical values of the parameters in the coupling condi- 
tions are shown in Table 1. A Rossby-Haurwitz wave with n = m = 1, 
A = 3 x 1( ) 7 m 2 s~ l , and U = 20 ms -1 , together with standard physical con- 
stants, was used to derive the data in this table. From these values it can be 
seen that Eq. 32 is satisfied, but Eq. 33 is certainly violated on intermediate 
and coarse grids. In terms of constructing a good smoother for the system, 
this means that Eq. 22 can be relaxed as though it were decoupled from 
the system, but the two remaining equations must be dealt with together, at 
least on coarse grids. In practise, it is easiest to use the same smoother on 
all grids to start with. 

To deal with Eqs. 22 and 23 together, collective relaxation is used. For 
linear equations, this means that, when the equations are relaxed at a point, 



corrections are made to all the variables associated with the equations such 
that the residuals of the equations become zero at that point. This may be 
done by replacing <f> and \ with (j) + 6$ and x + respectively, in Eqs. 23 
and 24 at a single point and the differential operators with their discretized 
counterparts and solving for the corrections 6$ and 8 X . Because Eq. 24 is 
nonlinear, a term proportional to 6<p6 x appears. We neglect this term and 
solve the resulting linear system directly. This method is equivalent to taking 
a single Newton step for these equations. 


4 Preliminary Numerical Results 

A preliminary code has been implemented that applies the multigrid method 
just described to the static equations. Eq. 22 was relaxed by red-black 
Gauss-Seidel iteration, and Eqs. 23 and 24 were relaxed collectively as de- 
scribed above in a lexicographic ordering. The equations are nonlinear, so 
the Full Approximation Scheme (FAS) [9] was used for the coarse- to- fine cor- 
rections. Full weighting was used for the fine to coarse grid restrictions, and 
linear interpolation for the coarse to fine grid transfers. Note that the grid 
transfers are straightforward because all the variables are defined on the same 
vertex centered grid. The standard five-point discretization was used for the 
Laplacian operator. Similarly, other derivatives were discretized using the 
usual finite difference formulae. At the time of writing, a semi-Lagrangian 
time-stepping scheme had been implemented, but the two codes had not been 
fully combined. 

In order to test the convergence of the multigrid scheme, we set the forc- 
ing functions in the static equations to a variety of functional forms. The 
magnitude of these functions was indicated by the Rossby-Haurwitz wave 
solution introduced in the previous section. When the problem was solved 
on a 64 x 32 grid with a V( 1 , 1 ) cycle, the convergence rates for the L 2 norm 
of the residuals were 0.22, 0.25, and 0.27 for Eqs. 22-24, respectively. When 
a V(2,l) cycle was used , the rates were 0.15, 0.13, and 0.14. In each of these 
cases, a single relaxation sweep consisted of relaxing Eq. 22 once followed by 
relaxing Eqs. 23 and 24 collectively once. 

These results suggest that multigrid may be an efficient way of solving 
these equations. Clearly, there are many possible variants on the scheme 
described above. For instance, the coupling analysis suggests that it may be 



fruitful to relax Eqs. 23 and 24 independently on fine grids and switch to 
collective relaxation only on coarser grids. 
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